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1. Introduction 

QCD at vanishing chemical potential and at increasing temperatures (T) undergoes a tran- 
sition at T=T C . In the cold, hadronic phase the dominant degrees of freedom are colourless 
hadrons, whereas in the high temperature, plasma phase the dominant degrees of free- 
dom are coloured objects (though the type of the transition -first order, second order or 
crossover- is not unambiguously determined yet, we use the expression 'phase' to charac- 
terize the system with different dominant degrees of freedom) . The equilibrium description 
as a function of the temperature is given by the equation of state (EoS). The complete 
determination of the EoS needs non-perturbative inputs, out of which lattice simulations 
is the most systematic approach. In lattice simulations there are two obvious difficulties, 
which are connected with 1.) the physical quark-mass limit and 2.) the continuum limit. 

ad 1.) The masses of the up and down quarks (m uc i) are small compared to the QCD scale. 
Since the computational costs increase sharply for small masses, most of the simula- 
tions are done at a larger mass than the physical one. 

ad 2.) The discretization of space-time through the lattice introduces the lattice spacing 
'a'. The continuum results are obtained in the a — >0 limit. Again, the computational 
costs increase sharply for small lattice spacings (for the EoS it scales approximately 
as 1/a 13 , which is a stronger growth of CPU-costs, than for typical T=0 simulations 
e.g. spectroscopy). In order to approach the continuum limit fast, improved actions 
are used. They all give the same QCD action in the continuum limit. However, 
for non-vanishing lattice spacing they have much less discretization artefacts than 
the most straightforward unimproved actions. In this paper we will use a Symanzik 
improved gauge and stout link improved fermionic action. 

The EoS has been determined in the continuum limit for the pure gauge theory (in this 
case -quenched simulations- the simulations are particularly easy as there is no fermionic 
degree of freedom. The simulated systems are equivalent to those where all the fermions 
are infinitely heavy, thus, is far from the physical situation) 1 [1, 2, 3]. Far less is known 

1 Note, however, that even in this relatively simple case there is still a few % difference between the 
different approaches. 
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about the unquenched case (QCD with dynamical quarks). There are published results 
for two-flavor QCD using unimproved staggered [4, 5], and improved Wilson fermions [6]. 
Few results are available for the 2+1 flavour case, among which the study done by Karsch, 
Laermann and Peikert in the year 2000 [7], using p4-improved staggered fermions, has often 
been used as the best result of EoS from lattice QCD. 

In this letter we calculate the EoS of dynamical QCD and attempt to improve on the 
result of Karsch, Laermann and Peikert by several means. 2 

1. ) We use for the lightest hadronic degree of freedom the physical pion with mass of 

mn «140 MeV instead of having the unphysical «600 MeV of Karsch, Laermann 
and Peikert in Ref. [7]. 

2. ) We use finer lattices than they [7], since our maximal temporal extension is Nt=6 

instead of their temporal extension of N t =4. 

3. ) As suggested for staggered QCD thermodynamics [8, 9], we keep our system on the 

line of constant physics (LCP) instead of increasing the physical quark mass with 
the temperature (if Karsch, Laermann and Peikert had cooled down e.g. two of their 
systems one at T=3T C and one at T=0.7T C down to T=0, the first system would 
have had approximately 4 times larger quark masses -two times larger pion masses- 
than the second one; this unphysical choice is known to lead systematics, which are 
comparable to the difference between the interacting and non-interacting plasma). 

4. ) For the staggered formulation of quarks the physically almost degenerate pion triplet 

has an unphysical non-degeneracy (so-called taste violation). This mass splitting 
Am 2 vanishes in the continuum limit as a — >0. Due to our smaller lattice spacing 
and particularly due to our stout-link improved [10] action the splitting Am 2 is much 
smaller than their splitting in Ref. [7]. 

5. ) In their Ref. [7] they used the approximate R-algorithm [11] in their simulations. 

This algorithm has an intrinsic parameter, the stepsize which, similarly to the lattice 
spacing, has to be extrapolated to zero. None of the previous staggered lattice ther- 
modynamic studies carried out this extrapolation. Using the R-algorithm without 
stepsize extrapolation leads to uncontrolled systematic errors. Instead of using the 
approximate R-algorithm this work uses the exact RHMC-algorithm (rational hybrid 
Monte-Carlo) [12, 13]. 

6. ) They used [7] the string tension to set the physical scale. This quantity, strictly speak- 

ing, does not exist in dynamical QCD (the string breaks and two mesons are produced 
[14]). Instead of this somewhat problematic quantity we use the quark-antiquark po- 
tential at intermediate distances. (This part of the potentials is theoretically sound 
[15] and it has been studied in the continuum limit [16]). 

For completeness -and for fairness- one should also mention that the analysis of [7] 
has also advantages (though as it will be shown, these advantages are also reached by the 
present letter in some way or another). 

2 It is important to emphasize, that these improvements are partly due to the increased CPU resources 
and partly due to theoretical improvements over the last 5 years. 
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1. ) At infinitely large T both the action of [7] and that of the present work approaches the 

continuum result as l/Nj*, though the prefactor for the action used by [7] is smaller. 
Note, however, that an extrapolation based on our N t =4 and 6 temporal extensions 
would give even in this limit roughly the same deviation from the continuum value 
as the action of Karsch, Laermann and Peikert at their N t =4 temporal extension. 

2. ) The thermodynamic limit is usually approached by large volumes V, which can be 

characterized by the aspect ratio r=N s /Nt (N s is the spatial extension of the lattice). 
The aspect ratio of Ref. [7] is 4. This is larger than r of the present analysis, which 
was in most of the cases 3. Note; however, that at several points we used aspect 
ratios of 3 and 4 in order to check, that the uncertainties due to r=3 are smaller than 
our statistical errors. 

In addition there is -at least- one uncertainty that neither [7] nor this work can really 
address. When determining the EoS at fixed temporal extension (in [7] with A^=4 and 
in this letter with A^=4,6) the lattice spacing changes when we change the temperature. 
Thus, the finite lattice spacing effects are different at low and at high temperatures. These 
uncertainties can be only resolved by continuum extrapolation. Though we have the EoS 
on two different sets of lattice spacings and one might attempt to do this extrapolation, it is 
fair to say that another set of lattice spacings is needed (N t =8). One of the reasons is, that 
in the hadronic phase, where the integration for the pressure starts, the lattice spacing 
is larger than 0.3 fm. In this region the lattice artefacts can not be really controlled 
(and in this deeply hadronic case it does not really help that an action is very good at 
asymptotically high temperatures, in the free non-interacting gas limit [17], as the p4 
action of [7] is). 

There are other ongoing thermodynamics projects (c.f. [18, 19]) improving on previous 
results. Both the MILC and the RBC-Bielefeld collaborations presented new results at the 
lattice conference. The MILC collaboration studies the equation of state along LCP's 
with two light quark masses (0.1 and 0.2 times m s ) at N t =4 and 6 lattices using asqtad 
improved staggered fermionic action [18]. The RBC-Bielefeld Collaboration developed the 
p4fat7 action to reduce taste symmetry violation. They study the Nf = 3 phase transition 
with several quark masses down to a pion mass of f» 320 MeV on N t =4 and 6 lattices [19]. 
Both works use the R algorithm. 

The paper is organized as follows. In Section 2 our definition for the Symanzik im- 
proved gauge and for the stout-link improved fermionic action is presented. We show our 
result on the line of constant physics and discuss how it was obtained. We demonstrate 
that our choice of action has small taste violation, when compared to other actions. The 
advantages of our exact RHMC simulation algorithm are discussed, too. Some details on 
the simulation points are summarized. Section 3 shortly discusses the methods (c.f. [20]) 
to determine the EoS and presents results at our two lattice spacings for the pressure, 
energy density, entropy density, quark number susceptibilities and the speed of sound. In 
Section 4 we summarize and conclude. 

2. Lattice action, simulations and the line of constant physics 

This section contains several technical details. Readers, who are not interested in these 
details, can simply skip to the next section. First we give our definition for the Symanzik 
improved gauge and for the stout-link improved fermionic action. We demonstrate that 
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Figure 1: Pion mass splitting A T = (m! 2 — m 2 ^) /T 2 as a function of (m^/X^) 2 . The lattice spacings 
are the same as those at the finite temperature transition point. The mass of the Goldstone pion is 
denoted by m^, that of the first non-Goldstone mode is by m' w . The horizontal blue line corresponds 
to the physical value of (m w /T c ) 2 , where T c = 173MeV was assumed [21]. The taste violation of our 
stout-link improved action is much smaller than that of the unimproved action and even somewhat 
smaller than that of the asqtad action at the same N t . 



our choice of stout-link improved staggered fermionic action has small taste violation, when 
compared to other staggered actions used in the literature to determine the EoS of QCD. 
The advantages of our exact RHMC simulation algorithm are emphasized. We discuss the 
importance of the LCP and show how to determine it by simulating in the three-flavour 
theory and by using the pseudoscalar and vector meson masses. Some details on the 
simulation points are summarized. 

Isotropic lattice couplings are used, thus the lattice spacings are identical in all direc- 
tions. The lattice action we used has the following form: 

S = Sg + Sf, (2.1) 

s 9 = E § ( c ° E *0*) + c i E Kf^ (2.2) 

Sf = Y,^u d {x)[m St ° U %y + ^ud^y]~ l/2 VuM 

x,y 

+ rj s (x)WU StOU %y + "i.g^.d/)}, (2-3) 

where W** 1 , W^* 2 are real parts of the traces of the ordered products of link matrices 
along the lxl, 1x2 rectangles in the fx, v plane. The coefficients satisfy cq + 8c\ = 1 and 
ci = —1/12 for the tree- level Symanzik improved action. rj u d and r] s are the pseudofermion 
fields for u, d and s quarks. Ip[U stout ) is the four-flavor staggered Dirac matrix with stout- 
link improvement [10]. Let us also note here, that we use the 4th root trick in eq. (2.1), 
which might lead to problems of locality. 

Our staggered action at a given Nt yields the same limit for the pressure at infinite 
temperatures as the standard unimproved action. There are various techniques improving 
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Figure 2: The line of constant physics. The result was obtained by using the (f> and K masses 
(see text). The strange quark mass in lattice units is shown as a function of /3. In the rest of our 
analysis we use light quark masses of m u d=m s /25. 

the high temperature scaling. However, usually an extrapolation based on Nt and Nt + 2 
with standard staggered action gives a better high T behavior for the pressure than p4 
[17] or asqtad [22, 23] action with Nt. Since our choice of action is about an order of 
magnitude faster than e.g. p4, we decided to use this less impoved action, with which our 
CPU resources made it possible to study two lattice spacings (iVt=4 and 6). 

Staggered fermions have an unconvenient property: they violate taste symmetry at 
finite lattice spacing. Among other things this violation results in a splitting in the pion 
spectrum, which should vanish in the continuum limit. The stout-link improvement makes 
the staggered fermion taste symmetry violation small already at moderate lattice spacings. 
We found that a stout-smearing level of N smr =2 and smearing parameter of p=0.15 are the 
optimal values of the smearing procedure. In order to illustrate the advantage of the stout- 
link action Figure 1 compares the taste violation in different approaches of the literature, 
which were used to determine the EoS of QCD. Results on the pion mass splitting for 
unimproved (used by Ref. [4, 5]) 3 , p4 improved (used by Ref. [7, 24]), asqtad improved 
(used by [18, 25]) and stout-link improved (this work) staggered fermions are shown. The 
parameters were chosen to be the ones used by the different collaborations at the finite 
temperature transition point. 

In previous staggered analyses the gauge configurations were produced by the R- 
algorithm [11] at a given stepsize. These studies were carried out usually at one stepsize, 
which is 1/2 or 2/3 of the light quark mass. The stepsize is an intrinsic parameter of 
the algorithm, which has to be extrapolated to zero. None of the previous staggered lat- 
tice thermodynamics studies performed this extrapolation. Using the R-algorithm without 
stepsize extrapolation leads to uncontrolled systematic errors. E.g. let us look at the dif- 
ference (on Nt=6 lattices at intermediate (3) between the extrapolated plaquette value and 
the value obtained at stepsize which is 2/3 of the light quark mass. This difference is larger 
than the total contribution of the plaquette to the pressure. Clearly, such a technique can 
not be used. 

Instead of using the approximate R-algorithm this work uses the exact RHMC- 
algorithm (rational hybrid Monte-Carlo) [12, 13]. This technique approximates the frac- 

3 We performed simulations to obtain the MILC standard action value at m q = 0.0125, f3 c = 5.415. 
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tional powers of the Dirac operator by rational functions. Since the condition number 
of the Dirac operator changes as we change the mass, one should determine the optimal 
rational approximation for each quark mass. Note however, that this should be done only 
once, and the obtained parameters of these functions can be used in the entire configura- 
tion production. Our choices for the rational approximation were as good as few times the 
machine precision for the whole range of the eigenvalues of the Dirac operator. 

Let us discuss the determination of the LCP. The LCP is defined as relationships be- 
tween the bare lattice parameters {(3 and lattice bare quark masses m u d and m s ). These 
relationships express that the physics (e.g. mass ratios) remains constant, while chang- 
ing any of the parameters. It is important to emphasize that the LCP is unambiguous 
(independent of the physical quantities, which are used to define the above relationships) 
only in the continuum limit ((3 — ► oo). For our lattice spacings fixing some relationships 
to their physical values means that some other relationships will slightly deviate from the 
physical one. In thermodynamics the relevance of LCP comes into play when the temper- 
ature is changed by (3 parameter. Then adjusting the mass parameters (m u d and m s ) is 
an important issue, neglecting this in simulations can lead to several % error in the EoS. 

A particularly efficient (however only approximate, see later) way to obtain an LCP is 
by using simulations with three degenerate flavors with lattice quark mass m q . The leading 
order chiral perturbation theory implies the mass relation for ss mesons. The strange quark 
mass is tuned accordingly, as 

m 2 p S l m v\m q =m s = (2m 2 K - ml) /ml, (2.4) 

where m ps and my are the pseudoscalar and vector meson masses in the simulations with 
three degenerate quarks. The light quark mass is calculated using the ratio m u d = m s /25 
obtained by experimental mass input in the chiral perturbation theory [26]. We obtain 
m s {(3) as shown in Figure 2. 

Our approach using eq. (2.4) is appropriate if in the nj=2+l theory the vector meson 
mass depends only weakly on the light quark masses and the chiral perturbation theory 
for meson masses works upto the strange quark mass. After applying the LCP we cross- 
checked the obtained spectrum of the n/=2+l simulations. The resulting pseudoscalar 
(pion, kaon) and phi mass ratios agree with the experimental values on the 5-10% level. 
In dynamical simulations approximately 10% effect on vq is observed by changing the light 
quark masses from the strange mass to the physical limit. Thus, this slight mismatch of 
the spectrum results in a subdominant error on our overall scale, which is less or around 
2%. 

The determination of the EoS needs quite a few simulation points. Results are needed 
on finite temperature lattices (7V t =4 or 6) and on zero temperature lattices (iV$ 3> 4 or 6) 
at several [3 values (we used 16 different (3 values for Nt=4 and 14 values for A^=6). Since 
our goal is to determine the EoS for physical quark masses we have to determine quantities 
in this small physical quark mass limit (we call these (3 dependent bare light quark masses 
m ud (phys)). 

For our finite temperature simulations (A^=4,6) we used physical quark masses. The 
spatial sizes were always at least 3 times the temporal sizes. For the whole (3 range on 
N t = 4 we checked that by increasing the N s /N t ratio from 3 to 4 the results remained the 
same within our statistical uncertainties. 

In the chirally broken phase (our zero temperature simulations, thus lattices for which 
N t S> 4 or 6, belong always to this class) chiral perturbation theory can be used to ex- 
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Table 1: Summary of our simulation points. For the physical light quark masses (we call them 
mud(phys)) 25 times smaller values were taken than for the strange mass. T^0 simulations were 
performed with the above m s and (3 pairs, and at 5 different m u d values: {l,3,5,7,9}-m u d(phys). 
T=0 simulations were performed with the above m s and f3 pairs, but at 4 different m uc [ val- 
ues: {3,5,7,9}-m u( i(phys). The total number of trajectories divided by 100 arc collected in the # 
columns. The left column shows the N t =4, whereas the right column shows the N t =6 data. (For 
an explanation of our labeling see the text.) 

trapolate by a controlled manner to the physical light quark masses. Therefore for most of 
our simulation points 4 we used four pion masses (m^ «235, 300, 355 and 405 MeV), which 
were somewhat larger than the physical one. (To simplify our notation in the rest of this 
paper we label these points as 3,5,7 and 9 times m u d(phys).) It turns out that the chiral 
condensates at all the four points can be fitted by linear function of pion mass squared with 
good x 2 - (Later we will show, that only the chiral condensate is to be extrapolated to get 
the EoS at the physical quark mass.) The volumes were chosen in a way, that for three out 
of these four quark masses the spatial extentions of the lattices were approximately equal 
or larger than four times the correlation lengths of the pion channel. We checked for a few 
(3 values that increasing the spatial and/or temporal extensions of the lattices results in 
the same expectation values within our statistical uncertainties. (For 3-m u d{phys) values 
the spatial lengths of the lattices were only three times the correlation length of the pion 
channel. However, excluding this point from the extrapolations, the results do not change.) 

A detailed list of our simulation points at zero and at non-zero temperature lattices 
are summarized in Table 1. 

3. Equation of state 

In this section results for two sets of lattice spacings (iVi=4,6) for the pressure, energy 
density, entropy density, quark number susceptibilities and for the speed of sound are 
presented. 

4 In the ft — 3.0. .3. 4 range the T = simulations were carried out at m u d(phys). 
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We shortly review the integral technique to obtain the pressure [20]. For large homo- 
geneous systems the pressure is proportional to the logarithm of the partition function: 



pa 



rri -t 

^3 log Z(T, V) = ^3 log Z(N S , N t ;0, m q ). 



V/a 



(3-1) 



(Index 'q' refers to the ud and s flavors.) The volume and temperature are connected to 
the spatial and temporal extensions of the lattice: 



V = (N s a) 3 , 



T = 



N t a 



(3.2) 



The divergent zero-point energy has to be removed by subtracting the zero temperature 
(Nf — > oo) part of eq. (3.1). In practice the zero temperature subtraction is performed by 
using lattices with finite, but large N t (called 7V t0 , see Table 1). So the normalized pressure 
becomes: 



rp4 t 



^3 log Z(N 3 , N t ; 0, m q ) - log Z(N s0 , N t0 ; /3, m q ) 



(3.3) 



With usual Monte-Carlo techniques one cannot measure log Z directly, but only its deriva- 
tives with respect to the bare parameters of the lattice action. Having determined the 
partial derivatives one integrates in the multi-dimensional parameter space: 



V_ 

rp4 



^ = Nf 



(13, m q ) 



(Po,m q0 ) 



d((3, m q 



1 / 3 log Z/d/3 
N t N3 \d\og Z/dm q 



1 



diogz /dp 

d log Zq I dm q 



, (3.4) 



where Z/Z$ are shorthand notations for Z(N S , Nt)/Z(N S Q, Nto). Since the integrand is 
a gradient, the result is by definition independent of the integration path. We need the 
pressure along the LCP, thus it is convenient to measure the derivatives of log Z along the 
LCP and perform the integration over this line in the (3, m u d and m s parameter space. 
The lower limits of the integrations (indicated by (3q and m q o) were set sufficiently below 
the transition point. By this choice the pressure gets independent of the starting point (in 
other words it vanishes at small temperatures). In the case of 2 + 1 flavor staggered QCD 
the derivatives of log Z with respect to (3 and m q are proportional to the expectation value 
of the gauge action ((S g ) c.f. eq. (2.1)) and to the chiral condensates {{tpip q ))^ respectively. 
Eq. (3.4) can be rewritten appropriately and the pressure is given by (in this formula we 
write out explicitely the flavours): 



p(l3,m ud ,m a ) 
J (/3o,m ud o,m s o) 



d(P,m ud ,m s 



1 



N t Ni 



(-Sg//?) 

(Mud) 

(# s ) 



<-_S g //0>o 
(^Vw)o 
(# s )o 



, (3-5) 

where (• • • )o means averaging on a _/V^ • N t Q lattice. 

The integral method was originally introduced for the pure gauge case, for which 
the integral is one dimensional, it is performed along the (3 axis. Previous studies for 
staggered dynamical QCD (e.g. [5, 27, 7]) used a one-dimensional parameter space instead 
of performing it along the LCP. Note, that for full QCD the integration should be performed 
along a LCP path in a multi-dimensional parameter space. 
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Table 2: Summary of the results for the 2+1 flavor pressure, speed of sound and 1 flavor quark 
number susceptibility in the non-interacting Stefan-Boltzmann limit. e/T 4 is 3 times, whereas s/T 3 
is 4 times the normalized value of the pressure (p/T 4 ) in the Stefan-Boltzmann limit. The first two 
lines gives the results for 7V t =4,6 and the third line contains the results in the continuum (in the 
thermodynamic limit). 

Using appropriate thermodynamical relations one can obtain any thermal properties 
of the system. For example the energy density (e), entropy density (s) and speed of sound 
(eg) can be derived as 

e = T(dp/dT)-p, s= ( e + p )T, C ^ = % (3 ' 6) 

To be able to do theses derivatives one has to know the temperature along the LCP. Since 
the temperature is connected to the lattice spacing as T = (A^a) -1 , we need a reliable 
estimate on a. The lattice spacings at different points of the LCP are determined by 
first matching the static potentials for different [3 values at an intermediate distance for 
m ud = {3, 5}m u d(phys) quark masses, then extrapolating the results to the physical quark 
mass. Relating these distances to physical observables (determining the overall scale in 
physical units) will be the topic of a subsequent publication. We show the results as a 
function of T/T c . The transition temperature (T c ) is defined by the inflection point of the 
isospin number susceptibility (%/, see later). 

To get the energy density the literature usually uses another quantity, namely e-3p, 
which can be also directly measured on the lattice. In our analysis it turned out to be more 
appropriate to calculate first the pressure directly from the raw lattice data (eq. (3.5)) and 
then determine the energy density and other quantities from the pressure (eq. (3.6)). The 
reasons for that can be summarized as follows. As we discussed we perform T^O simula- 
tions with physical quark masses, whereas the subtraction terms from T=0 simulations are 
extrapolated from larger quark masses. This sort of extrapolation is adequate for the chiral 
condensates, for which chiral perturbation techniques work well. Thus, one can choose an 
integration path for the T=0 part of the pressure, which moves along a LCP at some larger 
fnud (e.g. 9 times m U( i(phys)) and then at fixed (3 goes down to the physical quark mass. 
No comparable analogous technique is available for the combination e-3p. 

We have also calculated the pressure for the larger quark masses. Plotting it as a 
function of the temperature the differences between them are significant. As a function of 
T/T c these differences are smaller, but still remain statistically significant in the 1.2...2.0T C 
region. Note that statements on the mass dependence are only qualitative since such an 
analysis requires the careful matching of the scales at different quark masses. It is non- 
trivial and will be the topic of a forthcoming paper. 

Let us present the results. In order to show how the different quantities scale with 
the lattice spacing we show always N t =4,6 results on the same plot. In addition, in order 
to make the relationship with the continuum limit more transparent we multiply the raw 
lattice results at finite temporal extensions (AT t =4,6) with c con t/cN t , where the c values are 
the results in the free non- interacting plasma (Stefan-Boltzmann limit). These c values 
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Figure 3: a.) The left panel shows the pressure p, as a function of the temperature. Both 
N t —4: (red, upper curve) and N t =6 (blue, lower curve) data are obtained along the LCP. They are 
normalized by T 4 and scaled by c con t/cM t (see text and Table 2). In order to lead the eye lines 
connect the data points, b.) The right panel is the energy density (e), red (upper) and blue (lower) 
for iV t =4 and 6 respectively. This result was obtained directly from the pressure. 




1 1.5 2 2.5 3 

T/T c 

Figure 4: The interaction measure, the values are normalized by c co „t/c/v t of the energy density. 
The labeling is the same as for Figure 3. 

are summarized in Table 2 for the pressure, speed of sound, and for the quark number 
susceptibility at iVf=4,6 and in the continuum limit. By this multiplication the lattice 
thermodynamic quantities should approach the continuum Stefan-Boltzmann values for 
extreme large temperatures. 

Table 3 contains our most important numerical results. We tabulated the raw and 
normalized pressure values for both lattice spacings and for all of our simulation points. 
This data set and eq. (3.6) were used to obtain the following figures. Figure 3 shows 
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p 


T/T c 


p/T A (raw) 


p/T A (scaled) 


P 


T/T c 


p/T 4 (raw) 


p/T A (scaled) 


3.000 


0.90 


0.12(0.02) 


0.07(0.01) 


3.450 


0.80 


0.07(0.11) 


0.05(0.08) 


3.150 


0.95 


0.32(0.07) 


0.19(0.04) 


3.500 


0.87 


0.23(0.11) 


0.15(0.08) 


3.250 


0.98 


0.59(0.10) 


0.34(0.06) 


3.550 


0.96 


0.59(0.12) 


0.39(0.08) 


3.275 


0.99 


0.73(0.10) 


0.42(0.06) 


3.575 


1.02 


0.91(0.12) 


0.60(0.08) 


3.300 


1.01 


0.91(0.10) 


0.52(0.06) 


3.600 


1.07 


1.29(0.13) 


0.86(0.08) 


3.325 


1.04 


1.13(0.10) 


0.65(0.06) 


3.625 


1.14 


1.69(0.13) 


1.12(0.09) 


3.350 


1.06 


1.39 0.09 


0.79 0.05 

V / 


3.650 


1.20 


2.10(0.14) 


1.40 0.09 


3.400 


1.14 


2.04(0.10) 


1.16(0.06) 


3.675 


1.28 


2.51(0.14) 


1.66(0.10) 


3.450 


1.23 


2.79(0.10) 


1.59(0.06) 


3.700 


1.35 


2.88(0.15) 


1.91(0.10) 


3.500 


1.34 


3.56(0.11) 


2.04(0.07) 


3.750 


1.52 


3.50(0.15) 


2.32(0.10) 


3.550 


1.49 


4.32(0.12) 


2.47(0.07) 


3.800 


1.70 


3.99(0.16) 


2.65(0.11) 


3.600 


1.66 


4.96(0.12) 


2.83(0.07) 


3.850 


1.90 


4.36(0.16) 


2.89(0.11) 


3.650 


1.86 


5.46(0.12) 


3.12(0.07) 


3.930 


2.24 


4.82(0.17) 


3.19(0.11) 


3.700 


2.09 


5.84(0.12) 


3.34(0.07) 


4.000 


2.55 


5.14(0.17) 


3.41(0.11) 


3.850 


2.93 


6.57(0.15) 


3.75(0.09) 










4.000 


3.93 


6.97(0.16) 


3.98(0.09) 











Table 3: Numerical values of the pressure for all of our simulation points. The left column shows 
the -/V t =4, whereas the right column shows the N t =6 data. Both the raw values and the ones scaled 
by c cont /c Nt are given. 



the equation of state on N t =4,6 lattices. The pressure (left panel) and e (right panel) are 
presented as a function of the temperature. The Stefan-Boltzmann limit is also shown. 
On Figure 4 the interaction measure (i.e. (e — 3p)/T 4 ) is plotted, it is normalized by the 
corresponding c con t/cN t value of the energy density. Figure 5 shows the entropy density 
(left panel) and the speed of sound (right panel), which can be obtained by using the 
pressure and energy density data (c.f. sT=e+p and (? s =dp/de) of the previous Figure 3. 
Clearly, the uncertainties of the pressure and those of the energy density cumulate in the 
speed of sound, therefore it is less precisely determined. 

Light and strange quark number susceptibilities (xud and Xs) are defined via [28] 

, (3.7) 



T 2 Nf dn 2 . 



m 9 =o 



where [i u d and (i s are the light and strange quark chemical potentials (in lattice units). 
With the help of the quark number operators 

1 d 

Qq = 7^7- logdet($>+ m q ), 
the susceptibilities can be written as 



T 2 - N 3 I W^=° + \ fl/i, / 



The first term is usually referred as disconnected, the second as connected part. The 
connected part of the light quark number susceptibility is 2 times the susceptibility of the 
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Figure 6: The isospin susceptibility (left panel, normalized by T 2 ) and the connected part of the 
strangeness susceptibility (right panel, normalized by T 2 ). The labeling is the same as for Figure 
3. 

isospin number (%/). It is presented on the left panel of Figure 6. For our statistics and 
evaluation method the disconnected parts are all consistent with zero and their value is 
far smaller than those of the connected parts. The right panel of Figure 6 contains the 
connected part of the strange number susceptibility. 

4. Summary, conclusion 

In this letter we presented results of a large scale numerical lattice study on the thermo- 
dynamics of QCD. 



- 12 - 



We determined the equation of state. Our analysis attempted to improve on existing 
analyses by several means. We used for the lightest hadronic degree of freedom the physical 
pion mass. We used finer lattices with two different sets of lattice spacings (iV t =4,6). We 
kept our system on the line of constant physics (LCP) instead of changing the physics with 
the temperature. Due to our smaller lattice spacing and particularly due to our stout-link 
improved fermionic action the unphysical pion mass splitting was much smaller than in any 
previous analysis. We used an exact calculation algorithm instead of an approximate one. 
Our scale was determined by a theoretically sound quantity and not based on the string 
tension. 

We presented results for the pressure, energy density, entropy density, speed of sound 
and on the isospin and strangeness susceptibilities. 

Since the finite lattice spacing effects are quite different for different temperatures a 
reliable continuum estimate can only be given if results on an even finer lattice (Nt=8) 
were obtained. This sort of analysis would be a major step towards the final results for the 
equation of state. 
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